home *** CD-ROM | disk | FTP | other *** search
/ Complete Linux / Complete Linux.iso / docs / devel / lisp / akcl_lin.z / akcl_lin / mp / mpi-linux.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-03-08  |  38.0 KB  |  2,032 lines
  1.  
  2. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  3. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  4. /*~                                    ~*/
  5. /*~               OPERATIONS DE BASE (NOYAU)            ~*/
  6. /*~                                    ~*/
  7. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  8. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  9.  
  10. /* MODIFIED by M. Frigo (7 Mar 1993) to work with the akcl port
  11.  * of linux.
  12.  */
  13.  
  14. #ifdef linux
  15. # include "genpari.h"
  16.  
  17. /* +2147483648 */
  18.  
  19. ulong ABS_MOST_NEGS[3]={0x01ff0003, 0x01000003,1<<31};
  20. #define adecaler(x,tetpil,anavma) (RAVYZARC=(GEN)(x),((RAVYZARC>=(GEN)anavma)&&(RAVYZARC<(GEN)tetpil)))
  21.  
  22. GEN   mulsi(x,y)
  23.      long x;
  24.      GEN y;
  25. {
  26.   long s=signe(y),ly=lgef(y),i;
  27.   GEN z;
  28.   
  29.   if((!x)||(!s)) return gzero;
  30.   if(x<0) {s= -s;x= -x;}
  31.   z=cgeti(ly+1);hiremainder=0;
  32.   for(i=ly-1;i>=2;i--) z[i+1]=addmul(x,y[i]);
  33.   if(hiremainder) {z[2]=hiremainder;setlgef(z,ly+1);}
  34.   else {avma+=4;z[1]=z[0]-1;z++;setlgef(z,ly);}
  35.   setsigne(z,s);return z;
  36. }
  37.  
  38. GEN stoi(x)
  39.      long x;
  40. {
  41.   GEN y;
  42.   
  43.   if(!x) return gzero;
  44.   y=cgeti(3);
  45.   if(x>0) {y[1]=0x1000003;y[2]=x;}
  46.   else{y[1]=0xff000003;y[2]= -x;}
  47.   return y;
  48. }
  49.  
  50. GEN cgetg(x,y)
  51.      long x,y;
  52. {
  53.   unsigned long p1;
  54.   GEN z;
  55.   
  56.   p1=avma-(((unsigned short)x)<<2);if(p1<bot) err(errpile);
  57.   avma=p1;z=(GEN)p1;z[0]=0x10000+x+(y<<24);
  58.   return z;
  59. }
  60.  
  61. GEN cgeti(x)
  62.      long x;
  63. {
  64.   unsigned long p1;
  65.   GEN z;
  66.   
  67.   p1=avma-4*x;if(p1<bot) err(errpile);
  68.   avma=p1;z=(GEN)p1;z[0]=0x1010000+x;
  69.   return z;
  70. }
  71.  
  72. GEN cgetr(x)
  73.      long x;
  74. {
  75.   unsigned long p1;
  76.   GEN z;
  77.   
  78.   p1=avma-4*x;if(p1<bot) err(errpile);
  79.   avma=p1;z=(GEN)p1;z[0]=0x2010000+x;
  80.   return z;
  81. }
  82.  
  83. GEN icopy(x)
  84.      GEN x;
  85. {
  86.   GEN y;
  87.   long lx=lgef(x),i;
  88.   
  89.   y=cgeti(lx);
  90.   for(i=1;i<lx;i++) y[i]=x[i];
  91.   return y;
  92. }
  93.  
  94. GEN rcopy(x)
  95.      GEN x;
  96. {
  97.   GEN y;
  98.   long lx=lg(x),i;
  99.   
  100.   y=cgetr(lx);
  101.   for(i=1;i<lx;i++) y[i]=x[i];
  102.   return y;
  103. }
  104.  
  105.  
  106. GEN negi(x)
  107.      GEN x;
  108. {
  109.   long s=signe(x);
  110.   GEN y;
  111.   
  112.   if(!s) return gzero;
  113.   y=icopy(x);setsigne(y,-s);
  114.   return y;
  115. }
  116.  
  117. GEN negr(x)
  118.      GEN x;
  119. {
  120.   GEN y;
  121.   
  122.   y=rcopy(x);setsigne(y,-signe(x));
  123.   return y;
  124. }
  125.  
  126.  
  127. GEN absi(x)
  128.      GEN x;
  129. {
  130.   GEN y;
  131.   long s=signe(x);
  132.   
  133.   if(!s) return gzero;
  134.   y=icopy(x);setsigne(y,1);return y;
  135. }
  136.  
  137. GEN absr(x)
  138.      GEN x;
  139. {
  140.   GEN y;
  141.   long s=signe(x);
  142.   
  143.   y=rcopy(x);
  144.   if(s) setsigne(y,1);
  145.   return y;
  146. }
  147.  
  148. int expi(x)
  149.      GEN x;
  150. {
  151.   long lx=x[1]&0xffff;
  152.   
  153.   return lx==2 ? -8388608 : ((lx-2)<<5)-bfffo(x[2])-1;
  154. }
  155.  
  156. long itos(x)
  157.      GEN x;
  158. {
  159.   long s=signe(x),p2;
  160.   unsigned long p1;
  161.   
  162.   if(!s) return 0;
  163.   if(lgef(x)>3) err(affer2);
  164.   p1=x[2];if(p1>=0x80000000) err(affer2);
  165.   p2=(s>0)?p1:(-((long)p1));return p2;
  166. }
  167.  
  168. void mpaff(x,y)
  169.      GEN x,y;
  170. {
  171.   long tx=typ(x),ty=typ(y);
  172.   if(tx==1)
  173.     {
  174.       if(ty==1) affii(x,y);else affir(x,y);
  175.     }
  176.   else
  177.     {
  178.       if(ty==1) affri(x,y);else affrr(x,y);
  179.     }
  180. }
  181.  
  182. void affsi(s,x)
  183.      long s;
  184.      GEN x;
  185. {
  186.   long lx;
  187.   
  188.   if(!s) {x[1]=2;return;}
  189.   lx=lg(x);if(lx<3) err(affer1);
  190.   if(s>0) {x[1]=0x1000003;x[2]=s;}
  191.   else {x[1]=0xff000003;x[2]= -s;}
  192. }
  193.  
  194. void affsr(s,x)
  195.      long s;
  196.      GEN x;
  197. {
  198.   long l,i,d;
  199.   
  200.   if(!s)
  201.     {
  202.       l=((x[0]&0xffff)-2)<<5;x[1]=0x800000-l;x[2]=0;
  203.     }
  204.   else
  205.     {
  206.       d=1;if(s<0) {d= -1;s= -s;}
  207.       l=bfffo(s);setexpo(x,31-l);setsigne(x,d);
  208.       x[2]=(s<<l);for(i=3;i<lg(x);i++) x[i]=0;
  209.     }
  210. }
  211.  
  212. void affii(x,y)
  213.      GEN x,y;
  214. {
  215.   long lx=lgef(x),i;
  216.   
  217.   if(x==y) return;
  218.   if(lg(y)<lx) err(affer3);
  219.   for(i=1;i<lx;i++) y[i]=x[i];
  220. }
  221.  
  222. void affir(x,y)
  223.      GEN x,y;
  224. {
  225.   long lx=lgef(x),ly=lg(y),s,i,l,k;
  226.   
  227.   s=signe(x);
  228.   if(!s)
  229.     {
  230.       y[1]=0x800000-((ly-2)<<5);y[2]=0;
  231.     }
  232.   else
  233.     {
  234.       setsigne(y,s);l=((x[1]&0xffff)-2)<<5;s=bfffo(x[2]);
  235.       if(s)
  236.     {
  237.       setexpo(y,l-s-1);
  238.       if(lx<=ly)
  239.         {
  240.           for(i=lx;i<ly;i++) y[i]=0;k=0;
  241.           for(i=lx-1;i>=2;i--) {y[i]=shiftl(x[i],s)+k;k=hiremainder;}
  242.         }
  243.       else
  244.         {
  245.           shiftl(x[ly],s);k=hiremainder;
  246.           for(i=ly-1;i>=2;i--) {y[i]=shiftl(x[i],s)+k;k=hiremainder;}
  247.         }
  248.     }
  249.       else
  250.     {
  251.       setexpo(y,l-1);
  252.       if(lx<=ly)
  253.         {
  254.           for(i=lx;i<ly;i++) y[i]=0;
  255.           for(i=lx-1;i>=2;i--) y[i]=x[i];
  256.         }
  257.       else for(i=ly-1;i>=2;i--) y[i]=x[i];
  258.     }
  259.     }
  260. }
  261.  
  262. void affrr(x,y)
  263.      GEN x,y;
  264. {
  265.   long lx=lg(x),ly=lg(y),i;
  266.   
  267.   if(x==y) return;
  268.   y[1]=x[1];
  269.   if(!(x[1]&0xff000000)) y[2]=0;
  270.   else
  271.     {
  272.       if(lx<=ly)
  273.     {
  274.       for(i=2;i<lx;i++) y[i]=x[i];
  275.       for(i=lx;i<ly;i++) y[i]=0;
  276.     }
  277.       else for(i=2;i<ly;i++) y[i]=x[i];
  278.     }
  279. }
  280.  
  281. GEN shifts(x,y)
  282.      long x,y;
  283. {
  284.   long t[3];
  285.   
  286.   if(!x) return gzero;
  287.   t[0]=0x1010003;
  288.   if(x>0) {t[1]=0x1000003;t[2]=x;}
  289.   else {t[1]=0xff000003;t[2]= -x;}
  290.   return shifti(t,y);
  291. }
  292.  
  293. GEN shifti(x,n)
  294.      GEN x;
  295.      long n;
  296. {
  297.   long lx=lgef(x),i,s=signe(x),d,m,p1,p2,k;
  298.   GEN y;
  299.   
  300.   if(!s) return gzero;
  301.   if(!n) return icopy(x);
  302.   if(n>0)
  303.     {
  304.       d=n>>5;m=n&31;
  305.       if(m)
  306.     {
  307.       p1=shiftl(x[2],m);p2=hiremainder;k=0;
  308.       if(p2)
  309.         {
  310.           y=cgeti(lx+d+1);for(i=lx+1;i<=lx+d;i++) y[i]=0;
  311.           for(i=lx;i>=4;i--) {y[i]=shiftl(x[i-1],m)+k;k=hiremainder;}
  312.           y[3]=p1+k;y[2]=p2;
  313.         }
  314.       else
  315.         {
  316.           y=cgeti(lx+d);for(i=lx;i<lx+d;i++) y[i]=0;
  317.           for(i=lx-1;i>=3;i--) {y[i]=shiftl(x[i],m)+k;k=hiremainder;}
  318.           y[2]=p1+k;
  319.         }
  320.     }
  321.       else
  322.     {
  323.       y=cgeti(lx+d);for(i=lx;i<lx+d;i++) y[i]=0;
  324.       for(i=lx-1;i>=2;i--) y[i]=x[i];
  325.     }
  326.     }
  327.   else
  328.     {
  329.       n= -n;d=n>>5;m=n&31;if(lx<d+3) return gzero;
  330.       if(!m)
  331.     {
  332.       y=cgeti(lx-d);for(i=2;i<lx-d;i++) y[i]=x[i];
  333.     }
  334.       else 
  335.     {
  336.       m=32-m;d++;
  337.       p1=shiftl(x[2],m);
  338.       if(hiremainder)
  339.         {
  340.           y=cgeti(lx-d+1);y[2]=hiremainder;
  341.           for(i=3;i<=lx-d;i++)
  342.         {
  343.           p2=shiftl(x[i],m);y[i]=p1+hiremainder;p1=p2;
  344.         }
  345.         }
  346.       else
  347.         {
  348.           if(lx==d+2) return gzero;
  349.           y=cgeti(lx-d);
  350.           for(i=3;i<=lx-d;i++) 
  351.         {
  352.           p2=shiftl(x[i],m);y[i-1]=p1+hiremainder;p1=p2;
  353.         }
  354.         }
  355.     }
  356.     }   
  357.   y[1]=y[0];setsigne(y,s);return y;
  358. }
  359.  
  360. GEN shiftr(x,n)
  361.      GEN x;
  362.      long n;
  363. {
  364.   long l;
  365.   GEN y;
  366.   
  367.   y=rcopy(x);l=expo(x)+n;
  368.   if(l>0x7fffff||l<-0x800000) err(shier2);
  369.   setexpo(y,l);return y;
  370. }
  371.  
  372. GEN mptrunc(x)
  373.      GEN x;
  374. {
  375.   long e,i,s,lx=lg(x),p1,p2,m;
  376.   unsigned long d;
  377.   GEN y;
  378.   
  379.   if(typ(x)==1) return icopy(x);
  380.   s=signe(x);if(!s) return gzero;
  381.   e=expo(x);if(e<0) return gzero;
  382.   d=e>>5;m=e&31;if(d>=lx-2) err(truer2);
  383.   y=cgeti(d+3);y[1]=y[0];setsigne(y,s);
  384.   if(m==31) for(i=2;i<=d+2;i++) y[i]=x[i];
  385.   else
  386.     {
  387.       m++;p1=0;
  388.       for(i=2;i<=d+2;i++)
  389.     {
  390.       p2=shiftl(x[i],m);y[i]=hiremainder+p1;p1=p2;
  391.     }
  392.     }
  393.   return y;
  394. }
  395.  
  396. GEN mpent(x)
  397.      GEN x;
  398. {
  399.   long e,i,lx=lg(x),m,f,p1,p2;
  400.   unsigned long d;
  401.   GEN y,z;
  402.   
  403.   if(typ(x)==1) return icopy(x);
  404.   if(signe(x)>=0) return mptrunc(x);
  405.   e=expo(x);if(e<0) {y=cgeti(3);y[2]=1;y[1]=0xff000003;return y;}
  406.   d=e>>5;m=e&31;if(d>=lx-2) err(truer2);
  407.   y=cgeti(d+3);y[1]=0xff000003+d;
  408.   if(m==31) 
  409.     {
  410.       for(i=2;i<=d+2;i++) y[i]=x[i];
  411.       while((i<lx)&&(!x[i])) i++;
  412.       f=(i<lx);
  413.     }    
  414.   else
  415.     {
  416.       m++;p1=0;
  417.       for(i=2;i<=d+2;i++)
  418.     {
  419.       p2=shiftl(x[i],m);y[i]=hiremainder+p1;p1=p2;
  420.     }
  421.       if(p1) f=1;
  422.       else
  423.     {
  424.       while((i<lx)&&(!x[i])) i++;
  425.       f=(i<lx);
  426.     }
  427.     }
  428.   if(f)
  429.     {
  430.       for(i=d+2;(i>=2)&&(y[i]==0xffffffff);i--) y[i]=0;
  431.       if(i>=2) y[i]++;
  432.       else
  433.     {
  434.       z=y;y=cgeti(1);*y=(*z)+1;y[1]=z[1]+1;
  435.     }
  436.     }
  437.   return y;
  438. }
  439.  
  440. int mpcmp(x,y)
  441.      GEN x,y;
  442. {
  443.   if(typ(x)==1) return (typ(y)==1) ? cmpii(x,y) : cmpir(x,y);
  444.   return (typ(y)==1) ? -cmpir(y,x) : cmprr(x,y);
  445. }
  446.  
  447. int cmpsi(x,y)
  448.      long x;
  449.      GEN y;
  450. {
  451.   ulong p;
  452.   
  453.   if(!x) return -signe(y);
  454.   if(x>0)
  455.     {
  456.       if(signe(y)<=0) return 1;
  457.       if(lgef(y)>3) return -1;
  458.       p=y[2];if(p==x) return 0;
  459.       return (p<(ulong)x) ? 1 : -1;
  460.     }
  461.   else
  462.     {
  463.       if(signe(y)>=0) return -1;
  464.       if(lgef(y)>3) return 1;
  465.       p=y[2];if(p== -x) return 0;
  466.       return (p<(ulong)(-x)) ? -1 : 1;
  467.     }
  468. }
  469.  
  470. int cmpsr(x,y)
  471.      long x;
  472.      GEN y;      
  473. {
  474.   int p;
  475.   long av;
  476.   GEN z;
  477.   
  478.   if(!x) return -signe(y);
  479.   av=avma;z=cgetr(3);affsr(x,z);
  480.   p=cmprr(z,y);avma=av;return p;
  481. }    
  482.  
  483. int cmpii(x,y)
  484.      GEN x,y;
  485. {
  486.   long sx=signe(x),sy=signe(y),lx,ly,i;
  487.   
  488.   if(sx<sy) return -1;
  489.   if(sx>sy) return 1;
  490.   if(!sx) return 0;
  491.   lx=lgef(x);ly=lgef(y);
  492.   if(lx>ly) return sx;
  493.   if(lx<ly) return -sx;
  494.   for(i=2;(i<lx)&&(x[i]==y[i]);i++);
  495.   if(i==lx) return 0;
  496.   return ((ulong)x[i]>(ulong)y[i]) ? sx : -sx;
  497. }
  498.  
  499. int cmpir(x,y)
  500.      GEN x,y;
  501. {
  502.   long av=avma;
  503.   int p;
  504.   GEN z;
  505.   
  506.   if(!signe(x)) return -signe(y);
  507.   z=cgetr(lg(y));affir(x,z);
  508.   p=cmprr(z,y);avma=av;return p;
  509. }
  510.  
  511. int cmprr(x,y)
  512.      GEN x,y;
  513. {
  514.   long sx=signe(x),sy=signe(y),ex,ey,lx,ly,lz,i;
  515.   
  516.   if(sx<sy) return -1;
  517.   if(sx>sy) return 1;
  518.   if(!sx) return 0;
  519.   ex=expo(x);ey=expo(y);
  520.   if(ex>ey) return sx;
  521.   if(ex<ey) return -sx;
  522.   lx=lg(x);ly=lg(y);lz=(lx<ly)?lx:ly;
  523.   for(i=2;(i<lz)&&(x[i]==y[i]);i++);
  524.   if(i<lz) return ((ulong)x[i]>(ulong)y[i]) ? sx : -sx;
  525.   if(lx>=ly)
  526.     {
  527.       while((i<lx)&(!x[i])) i++;
  528.       return (i==lx) ? 0 : sx;
  529.     }
  530.   else
  531.     {
  532.       while((i<ly)&(!y[i])) i++;
  533.       return (i==ly) ? 0 : -sx;
  534.     }
  535. }    
  536.  
  537. GEN mpadd(x,y)
  538.      GEN x,y;
  539. {
  540.   if(typ(x)==1) return (typ(y)==1) ? addii(x,y) : addir(x,y);
  541.   return (typ(y)==1) ? addir(y,x) : addrr(x,y);
  542. }
  543.  
  544. GEN addss(x,y)
  545.      long x,y;
  546. {
  547.   long t[3];
  548.   
  549.   if(!x) return stoi(y);
  550.   t[0]=0x1010003;
  551.   if(x>0) {t[1]=0x1000003;t[2]=x;} else {t[1]=0xff000003;t[2]= -x;}
  552.   return addsi(y,t);
  553. }
  554.  
  555.  
  556. GEN addsi(x,y)
  557.      long x;
  558.      GEN y;
  559. {
  560.   long sx,sy,ly,p,i;
  561.   GEN z;
  562.   
  563.   if(!x) return icopy(y);
  564.   sy=signe(y);if(!sy) return stoi(x);
  565.   if(x<0) {sx= -1;x= -x;} else sx=1;
  566.   ly=lgef(y);
  567.   if(sx==sy)
  568.     {
  569.       p=addll(x,y[ly-1]);
  570.       if(overflow)
  571.     {
  572.       z=cgeti(ly+1);z[ly]=p;
  573.       for(i=ly-1;(i>2)&&(y[i-1]==0xffffffff);i--) z[i]=0;
  574.       if(i>2)
  575.         {
  576.           z[i]=y[i-1]+1;i--;while(i>=3) {z[i]=y[i-1];i--;}
  577.           z[2]=z[1]=z[0]-1;z++;avma+=4;
  578.         }
  579.       else {z[2]=1;z[1]=z[0];}
  580.     }
  581.       else
  582.     {
  583.       z=cgeti(ly);z[ly-1]=p;for(i=1;i<ly-1;i++) z[i]=y[i];
  584.     }
  585.       setsigne(z,sx);
  586.     }
  587.   else
  588.     {
  589.       if(ly==3)
  590.     {
  591.       if((ulong)y[2]>(ulong)x)
  592.         {
  593.           z=cgeti(3);z[1]=(sy<<24)+3;z[2]=y[2]-x;return z;
  594.         }
  595.       if(y[2]==x) return gzero;
  596.       z=cgeti(3);z[1]=((-sy)<<24)+3;z[2]=x-y[2];return z;
  597.     }
  598.       p=subll(y[ly-1],x);
  599.       if(overflow)
  600.     {
  601.       z=cgeti(ly);z[ly-1]=p;
  602.       for(i=ly-2;!(y[i]);i--) z[i]=0xffffffff;
  603.       z[i]=y[i]-1;
  604.       if((i>2)||z[i]) {i--;for(;i>=1;i--) z[i]=y[i];}
  605.       else
  606.         {
  607.           z[2]=z[1]=z[0]-1;z++;avma+=4;setsigne(z,sy);
  608.         }
  609.     }
  610.       else
  611.     {
  612.       z=cgeti(ly);z[ly-1]=p;for(i=1;i<ly-1;i++) z[i]=y[i];
  613.     }    
  614.     }
  615.   return z;
  616. }
  617.  
  618. GEN addii(x,y)
  619.      GEN x,y;
  620. {
  621.   long sx=signe(x),sy=signe(y),sz,lx=lgef(x),ly=lgef(y),i,j,p;
  622.   GEN z;
  623.   
  624.   if(!sx) return icopy(y);
  625.   if(!sy) return icopy(x);
  626.   if(lx<ly) {z=x;x=y;y=z;sz=sx;sx=sy;sy=sz;sz=lx;lx=ly;ly=sz;}
  627.   if(sx==sy)
  628.     {
  629.       z=cgeti(lx+1);overflow=0;
  630.       for(i=ly-1,j=lx-1;i>=2;i--,j--) z[j+1]=addllx(x[j],y[i]);
  631.       if(overflow)
  632.     {
  633.       for(;(j>=2)&&(x[j]==0xffffffff);j--) z[j+1]=0;
  634.       if(j>=2)
  635.         {
  636.           z[j+1]=x[j]+1;j--;
  637.           for(;j>=2;j--) z[j+1]=x[j];
  638.           z[1]=z[0]-1;z[2]=x[1];z++;avma+=4;
  639.         }
  640.       else {z[2]=1;z[1]=x[1]+1;}
  641.     }
  642.       else
  643.     {
  644.       for(;j>=2;j--) z[j+1]=x[j];
  645.       z[1]=z[0]-1;z[2]=x[1];z++;avma+=4;
  646.     }
  647.     }
  648.   else
  649.     {
  650.       if(lx==ly)
  651.     {
  652.       setsigne(y,1);setsigne(x,1);p=cmpii(x,y);
  653.       setsigne(y,sy);setsigne(x,sx);if(!p) return gzero;
  654.       if(p<0) {z=x;x=y;y=z;sz=sx;sx=sy;sy=sz;}
  655.     }
  656.       z=cgeti(lx);overflow=0;
  657.       for(i=ly-1,j=lx-1;i>=2;i--,j--) z[j]=subllx(x[j],y[i]);
  658.       if(overflow)
  659.     {
  660.       for(;!(x[j]);j--) z[j]=0xffffffff;
  661.       z[j]=x[j]-1;j--;
  662.       for(;j>=2;j--) z[j]=x[j];
  663.     }
  664.       else
  665.     {
  666.       for(;j>=2;j--) z[j]=x[j];
  667.     }
  668.       if(z[2]) z[1]=x[1];
  669.       else
  670.     {
  671.       for(j=3;(j<lx)&&(!z[j]);j++);
  672.       i=j-2;z[i+1]=z[i]=z[0]-i;z+=i;avma+=(i<<2);setsigne(z,sx);
  673.     }
  674.     }
  675.   return z;
  676. }      
  677.  
  678. GEN addsr(x,y)
  679.      long x;
  680.      GEN y;
  681. {
  682.   long p[3];
  683.   
  684.   if(!x) return rcopy(y);
  685.   p[0]=0x1010003;
  686.   if(x>0) {p[1]=0x1000003;p[2]=x;} else {p[1]=0xff000003;p[2]= -x;}
  687.   return addir(p,y);
  688. }
  689.  
  690. GEN addir(x,y)
  691.      GEN x,y;
  692. {
  693.   long l,e,ly,av,i,l1;
  694.   GEN z;
  695.   
  696.   if(!signe(x)) return rcopy(y);
  697.   if(!signe(y))
  698.     {
  699.       l=lgef(x)-(expo(y)>>5);if((l<3)||(l>32767)) err(adder3);
  700.       z=cgetr(l);affir(x,z);return z;
  701.     }
  702.   else
  703.     {
  704.       e=expo(y)-expi(x);ly=lg(y);
  705.       if(e>0)
  706.     {
  707.       l=ly-(e>>5);if(l<=2) return rcopy(y);
  708.     }
  709.       else
  710.     { 
  711.       l=ly+((-e)>>5)+1;if(l>32767) err(adder3);
  712.     }
  713.       av=avma;z=cgetr(l);affir(x,z);l1=av-avma;l=l1>>2;
  714.       z=addrr(z,y);
  715.       for(i=lg(z)-1;i>=0;i--) z[i+l]=z[i];z+=l;avma+=l1;
  716.     }
  717.   return z;
  718. }
  719.  
  720. GEN addrr(x,y)
  721.      GEN x,y;
  722. {
  723.   long sx=signe(x),sy=signe(y),lx=lg(x),ly=lg(y),lz,ex=expo(x),ey=expo(y),sz;
  724.   long av0=avma,e,l,i,d,m,flag,lp1,lp2,av,k,j,cex,f2;
  725.   GEN z,p1,p2;
  726.   
  727.   if(!sy)
  728.     {
  729.       if(!sx) {e=(ex>=ey)?ex:ey;z=cgetr(3);z[2]=0;z[1]=e+0x800000;return z;}
  730.       e=ex-ey;
  731.       if(e<=0) {z=cgetr(3);z[2]=0;z[1]=ey+0x800000;return z;}
  732.       l=(e>>5)+3;if(l>lx) l=lx;z=cgetr(l);
  733.       for(i=1;i<l;i++) z[i]=x[i];return z;
  734.     }
  735.   e=ey-ex;
  736.   if(!sx)
  737.     {
  738.       if(e<=0) {z=cgetr(3);z[2]=0;z[1]=ex+0x800000;return z;}
  739.       l=(e>>5)+3;if(l>ly) l=ly;z=cgetr(l);
  740.       for(i=1;i<l;i++) z[i]=y[i];return z;
  741.     }
  742.   if(e)
  743.     {
  744.       if(e<0) {z=x;x=y;y=z;lz=lx;lx=ly;ly=lz;ey=ex;e= -e;sz=sx;sx=sy;sy=sz;}
  745.       d=e>>5;m=e&31;
  746.       if(d>=ly-2) return rcopy(y);
  747.       l=d+lx;
  748.       if(l>=ly)
  749.     {
  750.       flag=1;p1=cgetr(ly);lp1=ly;lp2=ly-d;
  751.     }
  752.       else
  753.     {
  754.       flag=0;p1=cgetr(l+1);lp2=lx+1;lp1=l+1;
  755.     }
  756.       av=avma;
  757.       if(m)
  758.     {
  759.       p2=cgetr(lp2);m=32-m;
  760.       if(flag) {shiftl(x[lp2-1],m);k=hiremainder;}
  761.       else k=0;
  762.       for(i=lp2-1;i>=3;i--) 
  763.         {
  764.           p2[i]=shiftl(x[i-1],m)+k;k=hiremainder;
  765.         }
  766.       p2[2]=k;
  767.     }
  768.       else p2=x;
  769.     }
  770.   else
  771.     {
  772.       l=(lx>ly)?ly:lx;p1=cgetr(l);av=avma;lp2=lp1=l;flag=2;p2=x;m=0;
  773.     }
  774.   if(sx==sy)
  775.     {
  776.       overflow=0;
  777.       if(m+flag) for(i=lp1-1,j=lp2-1;j>=2;i--,j--) p1[i]=addllx(p2[j],y[i]);
  778.       else 
  779.     {
  780.       p1[lp1-1]=y[lp1-1];
  781.       for(i=lp1-2,j=lp2-2;j>=2;i--,j--) p1[i]=addllx(p2[j],y[i]);
  782.     }
  783.       if(overflow)
  784.     {
  785.       for(;(i>=2)&&(y[i]==0xffffffff);i--) p1[i]=0;
  786.       if(i>=2) {cex=0;p1[i]=y[i]+1;while(i>=3) {i--;p1[i]=y[i];}}
  787.       else 
  788.         {
  789.           cex=1;k=0x80000000;if(ey==0x7fffff) err(adder4);
  790.           for(i=2;i<lp1;i++) {p1[i]=shiftlr(p1[i],1)+k;k=hiremainder;}
  791.         }
  792.     }
  793.       else {cex=0;for(;i>=2;i--) p1[i]=y[i];}
  794.       p1[1]=(sx<<24)+ey+cex+0x800000;
  795.       avma=av;return p1;
  796.     }
  797.   else 
  798.     {
  799.       if(!e) 
  800.     {
  801.       for(i=2;(i<l)&&(p2[i]==y[i]);i++);
  802.       if(i==l)
  803.         {
  804.           e=ex-((l-2)<<5)+0x800000;if(e<0) err(adder5);
  805.           if(e>=0x1000000) err(adder4);
  806.           avma=av0;z=cgetr(3);z[2]=0;z[1]=e;return z;
  807.         }
  808.       else
  809.         {
  810.           f2=(((ulong)y[i])>((ulong)p2[i]))?1:0;
  811.         }
  812.     }
  813.       else f2=1;
  814.       if(f2)
  815.     {
  816.       overflow=0;
  817.       if(m+flag) for(i=lp1-1,j=lp2-1;j>=2;i--,j--) p1[i]=subllx(y[i],p2[j]);
  818.       else 
  819.         {
  820.           p1[lp1-1]=y[lp1-1];
  821.           for(i=lp1-2,j=lp2-2;j>=2;i--,j--) p1[i]=subllx(y[i],p2[j]);
  822.         }
  823.       if(overflow)
  824.         {
  825.           for(;(i>=2)&&(!y[i]);i--) p1[i]=0xffffffff;
  826.           p1[i]=y[i]-1;while(i>=3) {i--;p1[i]=y[i];}
  827.         }
  828.       else for(;i>=2;i--) p1[i]=y[i];
  829.     }
  830.       else
  831.     {
  832.       overflow=0;
  833.       if(m+flag) for(i=lp1-1;i>=2;i--) p1[i]=subllx(p2[i],y[i]);
  834.       else 
  835.         {
  836.           p1[lp1-1]=subllx(0,y[lp1-1]);
  837.           for(i=lp1-2;i>=2;i--) p1[i]=subllx(p2[i],y[i]);
  838.         }
  839.     }
  840.       for(i=2;!p1[i];i++);j=i-2;avma=av+(j<<2);p1[j]=p1[0]-j;p1+=j;
  841.       m=bfffo(p1[2]);e=ey-(j<<5)-m+0x800000;
  842.       if(e<0) err(adder5);
  843.       p1[1]=f2 ? (sy<<24)+e : (sx<<24)+e;
  844.       if(m)
  845.     {
  846.       k=0;for(i=lp1-1-j;i>=2;i--) {p1[i]=shiftl(p1[i],m)+k;k=hiremainder;}
  847.     }
  848.       return p1;
  849.     }
  850. }
  851.  
  852. GEN mpsub(x,y)
  853.      GEN x,y;
  854. {
  855.   if(typ(x)==1) return (typ(y)==1) ? subii(x,y) : subir(x,y);
  856.   return (typ(y)==1) ? subri(x,y) : subrr(x,y);
  857. }
  858.  
  859. GEN subii(x,y)
  860.      GEN x,y;
  861. {
  862.   long s=signe(y);
  863.   GEN z;
  864.   
  865.   if(x==y) return gzero;
  866.   setsigne(y,-s);z=addii(x,y);setsigne(y,s);
  867.   return z;
  868. }
  869.  
  870. GEN subrr(x,y)
  871.      GEN x,y;
  872. {
  873.   long s=signe(y);
  874.   GEN z;
  875.   
  876.   if(x==y)
  877.     {
  878.       z=cgetr(3);z[2]=0;z[1]=0x800000-(lg(x)<<5);return z;
  879.     }
  880.   setsigne(y,-s);z=addrr(x,y);setsigne(y,s);return z;
  881. }
  882.  
  883. GEN subsi(x,y)
  884.      long x;
  885.      GEN y;
  886. {
  887.   long s=signe(y);
  888.   GEN z;
  889.   
  890.   setsigne(y,-s);z=addsi(x,y);setsigne(y,s);return z;
  891. }
  892.  
  893. GEN subsr(x,y)
  894.      long x;
  895.      GEN y;
  896. {
  897.   long s=signe(y);
  898.   GEN z;
  899.   
  900.   setsigne(y,-s);z=addsr(x,y);setsigne(y,s);return z;
  901. }
  902.  
  903. GEN subss(x,y)
  904.      long x,y;
  905. {
  906.   return addss(-y,x);
  907. }
  908.  
  909.  
  910. GEN subir(x,y)
  911.      GEN x,y;
  912. {
  913.   long s=signe(y);
  914.   GEN z;
  915.   
  916.   setsigne(y,-s);z=addir(x,y);setsigne(y,s);return z;
  917. }
  918.  
  919. GEN subri(x,y)
  920.      GEN x,y;
  921. {
  922.   long s=signe(y);
  923.   GEN z;
  924.   
  925.   setsigne(y,-s);z=addir(y,x);setsigne(y,s);return z;
  926. }
  927.  
  928. GEN mpmul(x,y)
  929.      GEN x,y;
  930. {
  931.   if(typ(x)==1) return (typ(y)==1) ? mulii(x,y) : mulir(x,y);
  932.   return (typ(y)==1) ? mulir(y,x) : mulrr(x,y);
  933. }
  934.  
  935. GEN mulss(x,y)
  936.      long x,y;
  937. {
  938.   long s,p1;
  939.   GEN z;
  940.   
  941.   if((!x)||(!y)) return gzero;
  942.   s=1;if(x<0) {s= -1;x= -x;} if(y<0) {s= -s;y= -y;}
  943.   p1=mulll(x,y);
  944.   if(hiremainder) {z=cgeti(4);z[2]=hiremainder;z[3]=p1;}
  945.   else {z=cgeti(3);z[2]=p1;}
  946.   z[1]=z[0];setsigne(z,s);return z;
  947. }
  948.  
  949. GEN mulsr(x,y)
  950.      long x;
  951.      GEN y;
  952. {
  953.   long lx,i,k,s,p1,p2,e;
  954.   GEN z;
  955.   
  956.   if(!x) return gzero;
  957.   s=signe(y);if(x<0) {s= -s;x= -x;}
  958.   if(!s)
  959.     {
  960.       p1=bfffo(x);e=y[1]+31-p1;if(e>=0x1000000) err(muler2);
  961.       z=cgetr(3);z[1]=e;z[2]=0;
  962.     }
  963.   else
  964.     {
  965.       if(x==1) {z=rcopy(y);setsigne(z,s);return z;}
  966.       lx=lg(y);z=cgetr(lx);setsigne(z,s);
  967.       p2=mulll(x,y[lx-1]);
  968.       for(i=lx-2;i>=2;i--) z[i+1]=addmul(x,y[i]);
  969.       z[2]=hiremainder;p1=bfffo(hiremainder);
  970.       if(p1)
  971.     {
  972.       shiftl(p2,p1);k=hiremainder;
  973.       for(i=lx-1;i>=2;i--)
  974.         {
  975.           z[i]=shiftl(z[i],p1)+k;k=hiremainder;
  976.         }
  977.     }
  978.       e=32-p1+expo(y);if(e>=0x800000) err(muler2);
  979.       setexpo(z,e);
  980.     }  
  981.   return z;
  982. }
  983.  
  984. GEN mulii(x,y)
  985.      GEN x,y;
  986. {
  987.   long i,j,lx=lgef(x),ly=lgef(y),sx,sy,lz,p1,p2;
  988.   GEN z;
  989.   
  990.   sx=signe(x);if(!sx) return gzero;
  991.   sy=signe(y);if(!sy) return gzero;
  992.   if(sy<0) sx= -sx;
  993.   if(lx>ly) {z=x;x=y;y=z;lz=lx;lx=ly;ly=lz;}
  994.   lz=lx+ly-2;if(lz>=0x10000) err(muler1);
  995.   z=cgeti(lz);z[1]=z[0];setsigne(z,sx);
  996.   p1=x[lx-1];hiremainder=0;
  997.   for(i=ly-1;i>=2;i--) z[lx+i-2]=addmul(p1,y[i]);
  998.   z[lx-1]=hiremainder;
  999.   for(j=lx-2;j>=2;j--)
  1000.     {
  1001.       p1=x[j];hiremainder=0;
  1002.       for(i=ly-1;i>=2;i--)
  1003.     {
  1004.       p2=addmul(p1,y[i]);
  1005.       z[i+j-1]=addll(p2,z[i+j-1]);hiremainder+=overflow;
  1006.     }
  1007.       z[j]=hiremainder;
  1008.     }
  1009.   if(!(z[2]))
  1010.     {
  1011.       z[2]=z[1]-1;z[1]=z[0]-1;z++;avma+=4;
  1012.     }
  1013.   return z;
  1014. }
  1015.  
  1016. GEN mulrr(x,y)
  1017.      GEN x,y;
  1018. {
  1019.   long i,j,lx=lg(x),ly=lg(y),sx=signe(x),sy=signe(y),ex=expo(x),ey=expo(y);
  1020.   long e,flag,garde,p1,p2,lz;
  1021.   GEN z;
  1022.   
  1023.   e=ex+ey+0x800000;if(e>=0xffffff) err(muler4);
  1024.   if(e<0) err(muler5);
  1025.   if((!sx)||(!sy)) {z=cgetr(3);z[2]=0;z[1]=e;return z;}
  1026.   if(sy<0) sx= -sx;
  1027.   if(lx>ly) {lz=ly;z=x;x=y;y=z;flag=1;} else {lz=lx;flag=(lx!=ly);}
  1028.   z=cgetr(lz);z[1]=(sx<<24)+e;
  1029.   if(flag) mulll(x[2],y[lz]);else hiremainder=0;
  1030.   if(lz==3)
  1031.     {
  1032.       garde=flag ? addmul(x[2],y[2]) : mulll(x[2],y[2]);
  1033.       if((long)hiremainder<0) {z[2]=hiremainder;z[1]++;}
  1034.       else {z[2]=(garde<0)?(hiremainder<<1)+1:(hiremainder<<1);}
  1035.       return z;
  1036.     }
  1037.   p1=x[lz-1];garde=hiremainder;
  1038.   if(p1)
  1039.     {
  1040.       mulll(p1,y[3]);p2=addmul(p1,y[2]);
  1041.       garde=addll(p2,garde);z[lz-1]=overflow+hiremainder;
  1042.     }
  1043.   else z[lz-1]=0;
  1044.   for(j=lz-2;j>=3;j--)
  1045.     {
  1046.       p1=x[j];
  1047.       if(p1)
  1048.     {
  1049.       mulll(p1,y[lz+2-j]);
  1050.       p2=addmul(p1,y[lz+1-j]);
  1051.       garde=addll(p2,garde);hiremainder+=overflow;
  1052.       for(i=lz-j;i>=2;i--)
  1053.         {
  1054.           p2=addmul(p1,y[i]);
  1055.           z[i+j-1]=addll(p2,z[i+j-1]);hiremainder+=overflow;
  1056.         }
  1057.       z[j]=hiremainder;
  1058.     }
  1059.       else z[j]=0;
  1060.     }
  1061.   p1=x[2];p2=mulll(p1,y[lz-1]);
  1062.   garde=addll(p2,garde);hiremainder+=overflow;
  1063.   for(i=lz-2;i>=2;i--)
  1064.     {
  1065.       p2=addmul(p1,y[i]);
  1066.       z[i+1]=addll(p2,z[i+1]);hiremainder+=overflow;
  1067.     }
  1068.   z[2]=hiremainder;
  1069.   if((long)hiremainder>0)
  1070.     {
  1071.       overflow=(garde<0)?1:0;
  1072.       for(i=lz-1;i>=2;i--) {p1=z[i];z[i]=addllx(p1,p1);}
  1073.     }
  1074.   else z[1]++;
  1075.   return z;
  1076. }
  1077.  
  1078. GEN mulir(x,y)
  1079.      GEN x,y;
  1080. {
  1081.   long sx=signe(x),sy,av,lz,ey,e,garde,p1,p2,i,j;
  1082.   GEN z,temp;
  1083.   
  1084.   if(!sx) return gzero;
  1085.   sy=signe(y);ey=expo(y);
  1086.   if(!sy)
  1087.     {
  1088.       z=cgetr(3);z[2]=0;e=expi(x)+ey+0x800000;if(e>=0x1000000) err(muler6);
  1089.       z[1]=e;return z;
  1090.     }
  1091.   lz=lg(y);if(sy<0) sx= -sx;
  1092.   z=cgetr(lz);setsigne(z,sx);av=avma;
  1093.   temp=cgetr(lz+1);affir(x,temp);x=y;y=temp;
  1094.   e=expo(y)+ey+0x800000;if(e>=0xffffff) err(muler4);
  1095.   if(e<0) err(muler5);
  1096.   z[1]=(sx<<24)+e;
  1097.   mulll(x[2],y[lz]);
  1098.   if(lz==3)
  1099.     {
  1100.       garde=addmul(x[2],y[2]);
  1101.       if((long)hiremainder<0) {z[2]=hiremainder;z[1]++;}
  1102.       else {z[2]=(garde<0)?(hiremainder<<1)+1:(hiremainder<<1);}
  1103.       avma=av;return z;
  1104.     }
  1105.   garde=hiremainder;
  1106.   p1=x[lz-1];mulll(p1,y[3]);p2=addmul(p1,y[2]);
  1107.   garde=addll(p2,garde);z[lz-1]=overflow+hiremainder;
  1108.   for(j=lz-2;j>=3;j--)
  1109.     {
  1110.       p1=x[j];mulll(p1,y[lz+2-j]);
  1111.       p2=addmul(p1,y[lz+1-j]);
  1112.       garde=addll(p2,garde);hiremainder+=overflow;
  1113.       for(i=lz-j;i>=2;i--)
  1114.     {
  1115.       p2=addmul(p1,y[i]);
  1116.       z[i+j-1]=addll(p2,z[i+j-1]);hiremainder+=overflow;
  1117.     }
  1118.       z[j]=hiremainder;
  1119.     }
  1120.   p1=x[2];p2=mulll(p1,y[lz-1]);
  1121.   garde=addll(p2,garde);hiremainder+=overflow;
  1122.   for(i=lz-2;i>=2;i--)
  1123.     {
  1124.       p2=addmul(p1,y[i]);
  1125.       z[i+1]=addll(p2,z[i+1]);hiremainder+=overflow;
  1126.     }
  1127.   z[2]=hiremainder;
  1128.   if((long)hiremainder>0)
  1129.     {
  1130.       overflow=(garde<0)?1:0;
  1131.       for(i=lz-1;i>=2;i--) {p1=z[i];z[i]=addllx(p1,p1);}
  1132.     }
  1133.   else z[1]++;
  1134.   avma=av;return z;
  1135. }
  1136.  
  1137. GEN convi(x)
  1138.      GEN x;
  1139. {
  1140.   long lx,av=avma,lz;
  1141.   GEN z,p1,p2;  
  1142.   
  1143.   if(!signe(x))
  1144.     {
  1145.       z=cgeti(3);z[1]= -1;z[2]=0;avma=av;return z+3;
  1146.     }
  1147.   p1=absi(x);lx=lgef(p1);lz=((lx-2)*15)/14+3;z=cgeti(lz);z[1]= -1;
  1148.   for(p2=z+2;signe(p1);p2++) *p2=divisii(p1,1000000000,p1);
  1149.   avma=av;return p2;
  1150. }
  1151.  
  1152. GEN confrac(x)
  1153.      GEN x;
  1154. {
  1155.   long lx=lg(x),ex= -expo(x)-1,ex1,av=avma,ly,ey;
  1156.   long lr,nbdec,k,i,j;
  1157.   GEN y,res;
  1158.   
  1159.   ey=((lx-2)<<5)+ex;ly=(ey+63)>>5;y=cgeti(ly);ex1=ex>>5; /* 95 dans mp.s faux? */
  1160.   for(i=0;i<ex1;i++) y[i]=0;
  1161.   ex&=31;
  1162.   if(!ex) for(j=2;j<lx;j++) y[i++]=x[j];
  1163.   else
  1164.     {
  1165.       k=0;
  1166.       for(j=2;j<lx;j++) {y[i++]=shiftlr(x[j],ex)+k;k=hiremainder;}
  1167.       y[ly-2]=k;
  1168.     }
  1169.   y[ly-1]=0;
  1170.   nbdec=ey*0.30103+1;lr=(nbdec+17)/9;res=cgeti(lr);
  1171.   *res=nbdec;
  1172.   for(j=1;j<lr;j++)
  1173.     {
  1174.       hiremainder=0;
  1175.       for(i=ly-1;i>=0;i--) y[i]=addmul(y[i],1000000000);
  1176.       res[j]=hiremainder;
  1177.     }
  1178.   avma=av;return res;
  1179. }
  1180.  
  1181. void mulsii(x,y,z)
  1182.      long x;
  1183.      GEN y,z;
  1184. {
  1185.   long av=avma;
  1186.   GEN p1;
  1187.   
  1188.   p1=mulsi(x,y);affii(p1,z);avma=av;
  1189. }
  1190.  
  1191. void addsii(x,y,z)
  1192.      long x;
  1193.      GEN y,z;
  1194. {
  1195.   long av=avma;
  1196.   GEN p1;
  1197.   
  1198.   p1=addsi(x,y);affii(p1,z);avma=av;
  1199. }
  1200.  
  1201. long divisii(x,y,z)
  1202.      long y;
  1203.      GEN x,z;
  1204. {
  1205.   long av=avma,k;
  1206.   GEN p1;
  1207.   
  1208.   p1=divis(x,y);affii(p1,z);avma=av;
  1209.   k=hiremainder;return k;
  1210. }
  1211.  
  1212. long vals(x)
  1213.      long x;
  1214. {
  1215.   unsigned short int y,z;
  1216.   int s;
  1217.  
  1218.   if(!x) return -1;
  1219.   y=x;if(!y) {s=16;y=((ulong)x)>>16;} else s=0;
  1220.   z=y&255;if(!z) {s+=8;z=y>>8;}
  1221.   y=z&15;if(!y) {s+=4;y=z>>4;}
  1222.   z=y&3;if(!z) {s+=2;z=y>>2;}
  1223.   return (z&1) ? s : s+1;
  1224. }
  1225.  
  1226. long vali(x)
  1227.      GEN x;
  1228. {
  1229.   long i,lx=lgef(x);
  1230.   
  1231.   if(!signe(x)) return -1;
  1232.   for(i=lx-1;(i>=2)&&(!x[i]);i--);
  1233.   return ((lx-1-i)<<5)+vals(x[i]);
  1234. }
  1235.  
  1236. GEN mpdiv(x,y)
  1237.      GEN x,y;
  1238. {
  1239.   if(typ(x)==1) return (typ(y)==1) ? divii(x,y) : divir(x,y);
  1240.   return (typ(y)==1) ? divri(x,y) : divrr(x,y);
  1241. }
  1242.  
  1243. GEN divss(x,y)
  1244.      long x,y;
  1245. {
  1246.   long p1;
  1247.   
  1248.   if(!y) err(diver1);
  1249.   hiremainder=0;p1=divll((ulong)abs(x),(ulong)abs(y));
  1250.   if(y<0) {hiremainder= -((long)hiremainder);p1= -p1;}
  1251.   if(x<0) p1= -p1;
  1252.   return stoi(p1);
  1253. }
  1254.  
  1255. GEN modss(x,y)
  1256.      long x,y;
  1257. {
  1258.   long y1;
  1259.   
  1260.   if(!y) err(moder1);
  1261.   hiremainder=0;divll(abs(x),y1=abs(y));
  1262.   if(!hiremainder) return gzero;
  1263.   return (((long)hiremainder)<0) ? stoi(y1-hiremainder) : stoi(hiremainder);
  1264. }
  1265.  
  1266. GEN resss(x,y)
  1267.      long x,y;
  1268. {
  1269.   if(!y) err(reser1);
  1270.   hiremainder=0;divll(abs(x),abs(y));
  1271.   return (y<0) ? stoi(-((long)hiremainder)) : stoi(hiremainder);
  1272. }
  1273.  
  1274. GEN dvmdss(x,y,z)
  1275.      long x,y;
  1276.      GEN *z;
  1277. {
  1278.   GEN p1;
  1279.  
  1280.   p1=divss(x,y);*z=stoi(hiremainder);
  1281.   return p1;
  1282. }
  1283.  
  1284. void dvmdssz(x,y,z,t)
  1285.      long x,y;
  1286.      GEN z,t;
  1287. {
  1288.   long av=avma;
  1289.   GEN p1;
  1290.  
  1291.   p1=divss(x,y);affsi(hiremainder,t);
  1292.   mpaff(p1,z);avma=av;
  1293. }
  1294.  
  1295. GEN dvmdsi(x,y,z)
  1296.      long x;
  1297.      GEN y,*z;
  1298. {
  1299.   GEN p1;
  1300.   p1=divsi(x,y);*z=stoi(hiremainder);
  1301.   return p1;
  1302. }
  1303.  
  1304. void dvmdsiz(x,y,z,t)
  1305.      long x;
  1306.      GEN t,y,z;
  1307. {
  1308.   long av=avma;
  1309.   GEN p1;
  1310.   
  1311.   p1=divsi(x,y);affsi(hiremainder,t);
  1312.   mpaff(p1,z);avma=av;
  1313. }
  1314.  
  1315. GEN dvmdis(x,y,z)
  1316.      long y;
  1317.      GEN x,*z;
  1318. {
  1319.   GEN p1;
  1320.   p1=divis(x,y);*z=stoi(hiremainder);
  1321.   return p1;
  1322. }
  1323.  
  1324. void dvmdisz(x,y,z,t)
  1325.      long y;
  1326.      GEN x,t,z;
  1327. {
  1328.   long av=avma;
  1329.   GEN p1;
  1330.   
  1331.   p1=divis(x,y);affsi(hiremainder,t);
  1332.   mpaff(p1,z);avma=av;
  1333. }
  1334.  
  1335. void dvmdiiz(x,y,z,t)
  1336.      GEN x,y,z,t;
  1337. {
  1338.   long av=avma;
  1339.   GEN p1,p2;
  1340.  
  1341.   p1=dvmdii(x,y,&p2);mpaff(p1,z);mpaff(p2,t);
  1342.   avma=av;
  1343. }
  1344.  
  1345. GEN divsi(x,y)
  1346.      long x;
  1347.      GEN y;
  1348. {
  1349.   long s=signe(y),ly=lgef(y),p1;
  1350.  
  1351.   if(!s) err(diver2);
  1352.   if((!x)||(ly>3)||(y[2]<0)) {hiremainder=x;return gzero;}
  1353.   hiremainder=0;p1=divll(abs(x),y[2]);
  1354.   if(signe(y)<0) {hiremainder= -((long)hiremainder);p1= -p1;}
  1355.   if(x<0) p1= -p1;
  1356.   return stoi(p1);
  1357. }
  1358.  
  1359. GEN ressi(x,y)
  1360.      long x;
  1361.      GEN y;
  1362. {
  1363.   divsi(x,y);return stoi(hiremainder);
  1364. }
  1365.  
  1366. GEN modsi(x,y)
  1367.      long x;
  1368.      GEN y;
  1369. {
  1370.   long s;
  1371.   GEN p1;
  1372.   
  1373.   divsi(x,y);
  1374.   if(!hiremainder) return gzero;
  1375.   if(x>0) return stoi(hiremainder);
  1376.   else
  1377.     {
  1378.       s=signe(y);setsigne(y,1);p1=addsi(hiremainder,y);
  1379.       setsigne(y,s);return p1;
  1380.     }
  1381. }
  1382.  
  1383. GEN divis(y,x)
  1384.      long x;
  1385.      GEN y;
  1386. {
  1387.   long s=signe(y),ly=lgef(y),i,d;
  1388.   GEN z;
  1389.   
  1390.   if(!x) err(diver4);
  1391.   if(!s) {hiremainder=0;return gzero;}
  1392.   if(x<0) {s= -s;x= -x;} 
  1393.   if((ulong)x>(ulong)y[2])
  1394.     {
  1395.       if(ly==3) {hiremainder=itos(y);return gzero;}
  1396.       else {z=cgeti(ly-1);d=1;hiremainder=y[2];}
  1397.     }
  1398.   else {z=cgeti(ly);d=0;hiremainder=0;}
  1399.   for(i=d+2;i<ly;i++) z[i-d]=divll(y[i],x);
  1400.   z[1]=z[0];setsigne(z,s);if(s<0) hiremainder= -((long)hiremainder);
  1401.   return z;
  1402. }
  1403.  
  1404. GEN modis(x,y)
  1405.      long y;
  1406.      GEN x;
  1407. {
  1408.   divis(x,y);if(!hiremainder) return gzero;
  1409.   return (signe(x)>0) ? stoi(hiremainder) : stoi(abs(y)+hiremainder);
  1410. }
  1411.  
  1412. GEN resis(x,y)
  1413.      long y;
  1414.      GEN x;
  1415. {
  1416.   divis(x,y);return stoi(hiremainder);
  1417. }
  1418.      
  1419. GEN divir(x,y)
  1420.      GEN x,y;
  1421. {
  1422.   GEN xr,z;
  1423.   long av,ly;
  1424.   
  1425.   if(!signe(y)) err(diver5);
  1426.   if(!signe(x)) return gzero;
  1427.   ly=lg(y);z=cgetr(ly);av=avma;affir(x,xr=cgetr(ly+1));
  1428.   xr=divrr(xr,y);affrr(xr,z);avma=av;return z;
  1429. }
  1430.  
  1431. GEN divri(x,y)
  1432.      GEN x,y;
  1433.  
  1434. {
  1435.   GEN yr,z;
  1436.   long av,lx,ex,s=signe(y);
  1437.  
  1438.   if(!s) err(diver8);
  1439.   if(!signe(x))
  1440.   {
  1441.     ex=expo(x)-expi(y)+0x800000;
  1442.     if(ex<0) err(diver12);
  1443.     z=cgetr(3);z[1]=ex;z[2]=0;return z;
  1444.   }
  1445.   if((lg(y)==3)&&(y[2]>0)) return (s>0) ? divrs(x,y[2]) : divrs(x,-y[2]);
  1446.   lx=lg(x);z=cgetr(lx);av=avma;affir(y,yr=cgetr(lx+1));
  1447.   yr=divrr(x,yr);affrr(yr,z);avma=av;return z;
  1448. }
  1449.  
  1450. void diviiz(x,y,z)
  1451.      GEN x,y,z;
  1452. {
  1453.   long av=avma,lz;
  1454.   GEN p1,p2;
  1455.   
  1456.   if(typ(z)==1) {p1=divii(x,y);affii(p1,z);avma=av;}
  1457.   else
  1458.     {
  1459.       lz=lg(z);p1=cgetr(lz);p2=cgetr(lz);affir(x,p1);affir(y,p2);
  1460.       p1=divrr(p1,p2);affrr(p1,z);avma=av;
  1461.     }
  1462. }
  1463.      
  1464. void divrrz(x,y,z)
  1465.      GEN x,y,z;
  1466. {
  1467.   long av=avma;
  1468.   GEN p1;
  1469.  
  1470.   p1=divrr(x,y);mpaff(p1,z);avma=av;
  1471. }
  1472.  
  1473. void mpdivz(x,y,z)
  1474.      GEN x,y,z;
  1475. {
  1476.   long av=avma,lz;
  1477.   GEN p1,p2;
  1478.  
  1479.   if(typ(z)==1)
  1480.     {
  1481.       if(typ(x)==2||typ(y)==2) err(divzer1);
  1482.       p1=divii(x,y);affii(p1,z);avma=av;
  1483.     }
  1484.   else
  1485.     {
  1486.       if(typ(x)==1)
  1487.     {
  1488.       if(typ(y)==2) {p1=divir(x,y);mpaff(p1,z);avma=av;}
  1489.       else
  1490.         {
  1491.           lz=lg(z);p1=cgetr(lz);p2=cgetr(lz);affir(x,p1);affir(y,p2);
  1492.           p1=divrr(p1,p2);affrr(p1,z);avma=av;
  1493.         }
  1494.     }
  1495.       else
  1496.     {
  1497.       if(typ(y)==2) {p1=divrr(x,y);affrr(p1,z);avma=av;}
  1498.       else {p1=divri(x,y);affrr(p1,z);avma=av;}
  1499.     }
  1500.     }
  1501. }
  1502.  
  1503. GEN divsr(x,y)
  1504.      long x;
  1505.      GEN y;
  1506. {
  1507.   long av,ly;
  1508.   GEN p1,z;
  1509.  
  1510.   if(!signe(y)) err(diver3);
  1511.   if(!x) return gzero;
  1512.   ly=lg(y);z=cgetr(ly);av=avma;p1=cgetr(ly+1);affsr(x,p1);p1=divrr(p1,y);
  1513.   affrr(p1,z);avma=av;return z;
  1514. }
  1515.  
  1516. GEN modii(x,y)
  1517.      GEN x,y;
  1518. {
  1519.   long av=avma,tetpil;
  1520.   GEN p1;
  1521.  
  1522.   p1=dvmdii(x,y,-1);
  1523.   if(signe(p1)>=0) return p1;
  1524.   tetpil=avma;p1=(signe(y)>0) ? addii(p1,y) : subii(p1,y);
  1525.   return gerepile(av,tetpil,p1);
  1526. }
  1527.  
  1528. void modiiz(x,y,z)
  1529.      GEN x,y,z;
  1530. {
  1531.   long av=avma;
  1532.   GEN p1;
  1533.  
  1534.   p1=modii(x,y);affii(p1,z);avma=av;
  1535. }
  1536.  
  1537. void resiiz(x,y,z)
  1538.      GEN x,y,z;
  1539. {
  1540.   long av=avma;
  1541.   GEN p1;
  1542.  
  1543.   p1=resii(x,y);affii(p1,z);avma=av;
  1544. }
  1545.  
  1546. GEN divrs(x,y)
  1547.      long y;
  1548.      GEN x;
  1549. {
  1550.   long i,k,lx,ex,garde,sh,s=signe(x);
  1551.   GEN z;
  1552.  
  1553.   if(!y) err(diver6);
  1554.   if(!s)
  1555.     {
  1556.       z=cgetr(3);z[2]=0;z[1]=x[1]-31+bfffo(y);
  1557.       if(z[1]<0) err(diver7);return z;
  1558.     }
  1559.   if(y<0) {s= -s;y= -y;}
  1560.   if(y==1) {z=rcopy(x);setsigne(z,s);return z;}
  1561.   z=cgetr(lx=lg(x));setsigne(z,s);hiremainder=0;
  1562.   for(i=2;i<lx;i++) z[i]=divll(x[i],y);
  1563.   garde=divll(0,y);sh=bfffo(z[2]);ex=expo(x)-sh;if((-ex)>0x800000) err(diver7);
  1564.   setexpo(z,ex);shiftl(garde,sh);k=hiremainder;
  1565.   for(i=lx-1;i>=2;i--) {z[i]=shiftl(z[i],sh)+k;k=hiremainder;}
  1566.   return z;
  1567. }
  1568.  
  1569. mpdivis(x,y,z)
  1570.      GEN x,y,z;
  1571. {
  1572.   long av=avma;
  1573.   GEN p1,p2;
  1574.  
  1575.   p1=dvmdii(x,y,&p2);
  1576.   if(signe(p2)) {avma=av;return 0;}
  1577.   affii(p1,z);avma=av;return 1;
  1578. }
  1579.  
  1580. divise(x,y)
  1581.      GEN x,y;
  1582. {
  1583.   long av=avma;
  1584.   GEN p1;
  1585.  
  1586.   p1=dvmdii(x,y,-1);avma=av;
  1587.   return signe(p1) ? 0 : 1;
  1588. }
  1589.  
  1590.  
  1591. GEN dvmdii(x,y,z)
  1592.      GEN x,y,*z;
  1593. {
  1594.   long av,av2,lx,ly,lz,i,j,dec,sh,k,k1,sx=signe(x),sy=signe(y);
  1595.   long saux,k3,k4,av1,flk4;
  1596.   ulong si,qp;
  1597.   GEN p1,p2,p3,p4;
  1598.   
  1599.   if(!sy) err(dvmer1);
  1600.   if(!sx)
  1601.     {
  1602.       if(((long)z==0xffffffff)||((long)z==0)) return gzero;
  1603.       *z=gzero;return gzero;
  1604.     }
  1605.   lx=lgef(x);ly=lgef(y);lz=lx-ly;
  1606.   if(lz<0)
  1607.     {
  1608.       if((long)z==0xffffffff) return icopy(x);
  1609.       if(z==0) return gzero;
  1610.       *z=icopy(x);return gzero;
  1611.     }
  1612.   av=avma;if(sx<0) sy= -sy;
  1613.   if(ly==3)
  1614.     {
  1615.       si=y[2];
  1616.       if(si>(ulong)x[2]) {p1=cgeti(lx-1);hiremainder=x[2];dec=1;}
  1617.       else {p1=cgeti(lx);hiremainder=0;dec=0;}
  1618.       for(i=2+dec;i<lx;i++) p1[i-dec]=divll(x[i],si);
  1619.       if((long)z==0xffffffff)
  1620.     {
  1621.       avma=av;if(!hiremainder) return gzero;
  1622.       p2=cgeti(3);p2[1]=(sx<<24)+3;p2[2]=hiremainder;return p2;
  1623.     }
  1624.       if(lx!=(dec+2)) {p1[1]=p1[0];setsigne(p1,sy);} else {avma=av;p1=gzero;}
  1625.       if(z==0) return p1;
  1626.       if(!hiremainder) *z=gzero;
  1627.       else {p2=cgeti(3);p2[1]=(sx<<24)+3;p2[2]=hiremainder;*z=p2;}
  1628.       return p1;
  1629.     }
  1630.   else
  1631.     {
  1632.       p1=cgeti(lx);
  1633.       sh=bfffo(y[2]);
  1634.       if(sh)
  1635.     {
  1636.       p2=cgeti(ly);k=shiftl(y[2],sh);
  1637.       for(i=3;i<ly;i++) 
  1638.         {
  1639.           k1=shiftl(y[i],sh);p2[i-1]=k+hiremainder;k=k1;
  1640.         }
  1641.       p2[ly-1]=k;k=0;
  1642.       for(i=2;i<lx;i++)
  1643.         {
  1644.           k1=shiftl(x[i],sh);p1[i-1]=k+hiremainder;k=k1;
  1645.         }
  1646.       p1[lx-1]=k;
  1647.     }
  1648.       else {p1[1]=0;for(j=2;j<lx;j++) p1[j]=x[j];p2=y;}
  1649.       si=p2[2];saux=p2[3];
  1650.       for(i=1;i<=lz+1;i++)
  1651.     {
  1652.       if(p1[i]==si) 
  1653.         {
  1654.           qp=0xffffffff;k=addll(si,p1[i+1]);
  1655.         }
  1656.       else
  1657.         {
  1658.           hiremainder=p1[i];qp=divll(p1[i+1],si);
  1659.           overflow=0;k=hiremainder;
  1660.         }
  1661.       if(!overflow)
  1662.         {
  1663.           k1=mulll(qp,saux);k3=subll(k1,p1[i+2]);k+=overflow;
  1664.           flk4=((ulong)hiremainder>(ulong)k);k4=subll(hiremainder,k);
  1665.           while(flk4) {qp--;k3=subll(k3,saux);k4-=overflow;flk4=((ulong)k4>(ulong)si);k4=subll(k4,si);}
  1666.         }
  1667.       hiremainder=0;
  1668.       for(j=ly-1;j>=2;j--)
  1669.         {
  1670.           k1=addmul(qp,p2[j]);
  1671.           p1[i+j-1]=subll(p1[i+j-1],k1);hiremainder+=overflow;
  1672.         }
  1673.       if((ulong)p1[i]<(ulong)hiremainder)
  1674.         {
  1675.           overflow=0;qp--;
  1676.           for(j=ly-1;j>=2;j--) p1[i+j-1]=addllx(p1[i+j-1],p2[j]);
  1677.         }
  1678.       p1[i]=qp;
  1679.     }
  1680.       av1=avma;
  1681.       if((long)z!=0xffffffff)
  1682.     {
  1683.       if(p1[1])
  1684.         {
  1685.           p3=cgeti(lz+3);for(j=2;j<=lz+2;j++) p3[j]=p1[j-1];
  1686.         }
  1687.       else
  1688.         {
  1689.           p3=cgeti(lz+2);
  1690.           if(!lz) sy=0;else for(j=2;j<=lz+1;j++) p3[j]=p1[j];
  1691.         }
  1692.       if(lg(p3)<3) p3[1]=2;else {p3[1]=p3[0];setsigne(p3,sy);}
  1693.     }
  1694.       if(z!=0)
  1695.     {
  1696.       for(j=lz+2;(j<lx)&&(!p1[j]);j++);
  1697.       if(j==lx) p4=gzero;
  1698.       else
  1699.         {
  1700.           p4=cgeti(lx-j+2);p4[1]=p4[0];
  1701.           if(!sh) for(i=2;j<lx;j++,i++) p4[i]=p1[j];
  1702.           else
  1703.         {
  1704.           hiremainder=0;k1=shiftlr(p1[j++],sh);k=hiremainder;
  1705.           if(k1) {p4[2]=k1;dec=1;} else {p4[1]=p4[0]-1;p4++;avma+=4;p4[1]=p4[0];dec=0;}
  1706.           for(i=2+dec;j<lx;j++,i++)
  1707.             {
  1708.               p4[i]=shiftlr(p1[j],sh)+k;k=hiremainder;
  1709.             }
  1710.         }
  1711.           setsigne(p4,sx);
  1712.         }
  1713.     }
  1714.       if((long)z==0xffffffff) return gerepile(av,av1,p4);
  1715.       if((long)z==0) return gerepile(av,av1,p3);
  1716.       av2=avma;dec=lpile(av,av1,0)>>2;*z=adecaler(p4,av1,av2)?p4+dec:p4;
  1717.       return adecaler(p3,av1,av2)?p3+dec:p3;
  1718.     }
  1719. }
  1720.  
  1721. GEN divrr(x,y)
  1722.      GEN x,y;
  1723. {
  1724.   long sx=signe(x),sy=signe(y),lx,ly,lz,ex,ex1,i,z0;
  1725.   long ldif,y0,y1,si,saux,qp,k,k1,k3,k4,j,flk4;
  1726.   GEN z;
  1727.   
  1728.   if(!sy) err(diver9);
  1729.   ex=expo(x)-expo(y)+0x800000;
  1730.   if(ex<=0) err(diver10);
  1731.   if(ex>=0x1000000) err(diver11);
  1732.   if(!sx)
  1733.     {
  1734.       z=cgetr(3);z[1]=ex;z[2]=0;return z;
  1735.     }
  1736.   lx=lg(x);ly=lg(y);lz=(lx<=ly)?lx:ly;
  1737.   z=cgetr(lz);if(sy<0) sx= -sx;
  1738.   ex1=(sx<<24)+ex;
  1739.   if(ly==3)
  1740.     {
  1741.       i=x[2];si=(lx>3)?x[3]:0;
  1742.       if((ulong)i<(ulong)y[2])
  1743.     {
  1744.       hiremainder=i;z[2]=divll(si,y[2]);
  1745.       z[1]=ex1-1;return z;
  1746.     }
  1747.       else
  1748.     {
  1749.       hiremainder=((ulong)i)>>1;
  1750.       z[2]=(i&1)?divll((((ulong)si)>>1)|(0x80000000),y[2]):divll(((ulong)si)>>1,y[2]);
  1751.       z[1]=ex1;return z;
  1752.     }
  1753.     }
  1754.   z0= *z;*z=0;
  1755.   for(i=2;i<=lz-1;i++) z[i-1]=x[i];
  1756.   z[lz-1]=(lx>lz) ? x[lz] : 0;
  1757.   ldif=ly-lz;if(!ldif) {y0=y[lz];y[lz]=0;}
  1758.   if(ldif<=1) {y1=y[lz+1];y[lz+1]=0;}
  1759.   si=y[2];saux=y[3];
  1760.   for(i=0;i<lz-1;i++)
  1761.     {
  1762.       if(z[i]==si) 
  1763.     {
  1764.       qp=0xffffffff;k=addll(si,z[i+1]);
  1765.     }
  1766.       else
  1767.     {
  1768.       hiremainder=z[i];qp=divll(z[i+1],si);
  1769.       overflow=0;k=hiremainder;
  1770.     }
  1771.       if(!overflow)
  1772.     {
  1773.       k1=mulll(qp,saux);k3=subll(k1,z[i+2]);k+=overflow;
  1774.       flk4=((ulong)hiremainder>(ulong)k);k4=subll(hiremainder,k);
  1775.       while(flk4) {qp--;k3=subll(k3,saux);k4-=overflow;flk4=((ulong)k4>(ulong)si);k4=subll(k4,si);}
  1776.     }
  1777.       mulll(qp,y[lz+1-i]);
  1778.       for(j=lz-i;j>=2;j--)
  1779.     {
  1780.       k1=addmul(qp,y[j]);
  1781.       z[i+j-1]=subll(z[i+j-1],k1);hiremainder+=overflow;
  1782.     }
  1783.       if((ulong)z[i]<(ulong)hiremainder)
  1784.     {
  1785.       overflow=0;qp--;
  1786.       for(j=lz-i;j>=2;j--) z[i+j-1]=addllx(z[i+j-1],y[j]);
  1787.     }
  1788.       z[i]=qp;
  1789.     }
  1790.   if(!ldif) y[lz]=y0;if(ldif<=1) y[lz+1]=y1;
  1791.   for(j=lz-1;j>=2;j--) z[j]=z[j-1];
  1792.   if(*z)
  1793.     {
  1794.       k=0x80000000;
  1795.       for(j=2;j<lz;j++) {z[j]=shiftlr(z[j],1)+k;k=hiremainder;}
  1796.     }
  1797.   else ex1--;
  1798.   z[1]=ex1;*z=z0;return z;
  1799. }
  1800.  
  1801.   
  1802. GEN gerepile(l,p,q)
  1803.     GEN l,p,q;
  1804.  
  1805. {
  1806.   long av,declg,tl;
  1807.   GEN ll,pp,l1,l2,l3;
  1808.  
  1809.   declg=(long)l-(long)p;if(declg<=0) return q;
  1810.   for(ll=l,pp=p;pp>(GEN)avma;) *--ll= *--pp;
  1811.   av=(long)ll;
  1812.   while((ll<l)||((ll==l)&&(long)q))
  1813.   {
  1814.     l2=ll+lontyp[tl=typ(ll)];
  1815.     if(tl==10) {l3=ll+lgef(ll);ll+=lg(ll);if(l3>ll) l3=l2;}
  1816.     else {ll+=lg(ll);l3=ll;} 
  1817.     for(;l2<l3;l2++) 
  1818. /*    for(;l2<ll;l2++) */
  1819.       {
  1820.     l1=(GEN)(*l2);
  1821.     if((l1<l)&&(l1>=(GEN)avma))
  1822.       {
  1823.         if(l1<p) *l2+=declg;
  1824.         else
  1825.           if(ll<=l) err(gerper);
  1826.       }
  1827.       }
  1828.   }
  1829.   if((!((long)q))||((q<p)&&(q>=(GEN)avma)))
  1830.   {
  1831.     avma=av;return q+(declg>>2);
  1832.   }
  1833.   else {avma=av;return q;}
  1834. }
  1835.  
  1836. void cgiv(x)
  1837.      GEN x;
  1838. {
  1839.   long p;
  1840.  
  1841.   if((p=pere(x))==255) return;
  1842.   if((x!=(GEN)avma)||(p>1)) {setpere(x,p-1);return;}
  1843.   do x+=lg(x);while(!pere(x));
  1844.   avma=(long)x;
  1845.   return;
  1846. }
  1847.  
  1848. ulong overflow,hiremainder;
  1849.  
  1850. int addll(x,y)
  1851.      int x,y;
  1852. {
  1853.   int z;
  1854.  
  1855.   z=x+y;
  1856.   if((x>=0)&&(y>=0)) {overflow=0;return z;}
  1857.   if((x<0)&&(y<0)) {overflow=1;return z;}
  1858.   overflow=(z>=0)?1:0; return z;
  1859. }
  1860.  
  1861. int addllx(x,y)
  1862.      int x,y;
  1863. {
  1864.   int z;
  1865.  
  1866.   z=x+y+overflow;
  1867.   if((x>=0)&&(y>=0)) {overflow=0;return z;}
  1868.   if((x<0)&&(y<0)) {overflow=1;return z;}
  1869.   overflow=(z>=0)?1:0; return z;
  1870. }
  1871.  
  1872. int subll(x,y)
  1873.      int x,y;
  1874. {
  1875.   int z;
  1876.  
  1877.   z=x-y;
  1878.   if((x>=0)&&(y<0)) {overflow=1;return z;}
  1879.   if((x<0)&&(y>=0)) {overflow=0;return z;}
  1880.   overflow=(z>=0)?0:1; return z;
  1881. }
  1882.  
  1883. int subllx(x,y)
  1884.      int x,y;
  1885. {
  1886.   int z;
  1887.  
  1888.   z=x-y-overflow;
  1889.   if((x>=0)&&(y<0)) {overflow=1;return z;}
  1890.   if((x<0)&&(y>=0)) {overflow=0;return z;}
  1891.   overflow=(z>=0)?0:1; return z;
  1892. }
  1893.  
  1894. int shiftl(x,y)
  1895.      ulong x,y;
  1896. {
  1897.   hiremainder=x>>(32-y);return (x<<y);
  1898. }
  1899.  
  1900. int shiftlr(x,y)
  1901.      ulong x,y;
  1902. {
  1903.   hiremainder=x<<(32-y);return (x>>y);
  1904. }
  1905.  
  1906. int bfffo(x)
  1907.      ulong x;
  1908. {
  1909.   int sc;
  1910.   static int tabshi[16]={4,3,2,2,1,1,1,1,0,0,0,0,0,0,0,0};
  1911.  
  1912.   if(x&(0xffff0000)) sc=0;else {sc=16;x<<=16;}
  1913.   if(!(x&(0xff000000))) {sc+=8;x<<=8;}
  1914.   if(x&(0xf0000000)) x>>=28;else {sc+=4;x>>=24;}
  1915.   sc+=tabshi[x];return sc;
  1916. }
  1917.  
  1918. int mulll(x,y)
  1919.      ulong x,y;
  1920. {
  1921.   ulong xlo,xhi,ylo,yhi;
  1922.   ulong z;
  1923.  
  1924.   xlo=x&65535;xhi=x>>16;ylo=y&65535;yhi=y>>16;
  1925.   z=addll(xlo*yhi,xhi*ylo);
  1926.   hiremainder=(overflow)?xhi*yhi+65536+(z>>16):xhi*yhi+(z>>16);
  1927.   z=addll(xlo*ylo,(z<<16));hiremainder+=overflow;
  1928.   return z;
  1929. }
  1930.  
  1931. int addmul(x,y)
  1932.      ulong x,y;
  1933. {
  1934.   ulong xlo,xhi,ylo,yhi;
  1935.   ulong z,z2;
  1936.  
  1937.   xlo=x&65535;xhi=x>>16;ylo=y&65535;yhi=y>>16;
  1938.   z=addll(xlo*yhi,xhi*ylo);
  1939.   z2=(overflow)?xhi*yhi+65536+(z>>16):xhi*yhi+(z>>16);
  1940.   z=addll(xlo*ylo,(z<<16));z2+=overflow;
  1941.   z=addll(z,hiremainder);hiremainder=z2+overflow;
  1942.   return z;
  1943. }
  1944.  
  1945. /* Ancienne version de divll, en generale plus lente sauf si le processeur flottant
  1946. est tres rapide
  1947.  
  1948. int divll(x,y)
  1949.      ulong x,y;
  1950. {
  1951.   ulong p,p1,p2;
  1952.   long p3,p4;
  1953.   double dbl;
  1954.  
  1955.   if((ulong)hiremainder>=y) err(divller1);
  1956.   p1=hiremainder;dbl=(C4*p1+x)/y;
  1957.   if(dbl < C4/2) p=(ulong)dbl;
  1958.   else
  1959.     {
  1960.       if(dbl>C4-1) p=0xffffffff;
  1961.       else{dbl-=(C4/2);p=(ulong)dbl;p|=0x80000000;}
  1962.     }
  1963.   p2=mulll(p,y);p3=x-p2;
  1964.   if(x<p2) hiremainder++;
  1965.   if((p1==(ulong)hiremainder)&&((ulong)p3<y))
  1966.     {hiremainder=p3;return p;}
  1967.   if(p1<(ulong)hiremainder) {hiremainder=p3+y;return p-1;}
  1968.   hiremainder=p3-y;return p+1;
  1969. }
  1970.  
  1971. */
  1972.  
  1973. int divll(x,y)
  1974.      ulong x,y;
  1975. {
  1976. #define HIBIT 0x80000000
  1977. #define HIMASK 0xffff0000
  1978. #define LOMASK 0xffff
  1979. #define HIWORD(a) (a >> 16)
  1980. /* si le compilateur est bugge, il faut mettre (a >> 16) & LOMASK) */
  1981. #define LOWORD(a) (a & LOMASK)
  1982. #define GLUE(hi, lo) ((hi << 16) + lo)
  1983. #define SPLIT(a, b, c) b = HIWORD(a); c = LOWORD(a)
  1984.  
  1985.     ulong v1, v2, u3, u4, q1, q2, aux, aux1, aux2;
  1986.     int k;
  1987.     
  1988.     for(k = 0; !(y & HIBIT); k++)
  1989.         {
  1990.             hiremainder <<= 1;
  1991.             if (x & HIBIT) hiremainder++;
  1992.             x <<= 1;
  1993.             y <<= 1;
  1994.         }
  1995.         
  1996.     SPLIT(y, v1, v2);
  1997.     SPLIT(x, u3, u4);
  1998.     
  1999.     q1 = hiremainder / v1; if (q1 & HIMASK) q1 = LOMASK;
  2000.     hiremainder -= q1 * v1;
  2001.     aux = v2 * q1;
  2002. again:
  2003.     SPLIT(aux, aux1, aux2);
  2004.     if (aux2 > u3) aux1++;
  2005.     if (aux1 > hiremainder) {q1--; hiremainder += v1; aux -= v2; goto again;}
  2006.     u3 -= aux2;
  2007.     hiremainder -= aux1;
  2008.     hiremainder <<= 16; hiremainder += u3 & LOMASK;
  2009.     
  2010.     q2 = hiremainder / v1; if (q2 & HIMASK) q2 = LOMASK;
  2011.     hiremainder -= q2 * v1;
  2012.     aux = v2 * q2;
  2013. again2:
  2014.     SPLIT(aux, aux1, aux2);
  2015.     if (aux2 > u4) aux1++;
  2016.     if (aux1 > hiremainder) {q2--; hiremainder += v1; aux -= v2; goto again2;}
  2017.     u4 -= aux2;
  2018.     hiremainder -= aux1;
  2019.     hiremainder <<= 16; hiremainder += u4 & LOMASK;
  2020.     hiremainder >>= k;
  2021.     return GLUE(q1, q2);
  2022. }
  2023.  
  2024. long mulmodll(a,b,c)
  2025.      ulong a,b,c;
  2026. {
  2027.   divll(mulll(a,b),c);
  2028.   return hiremainder;
  2029. }
  2030.  
  2031. #endif
  2032.